Introduction

This markdown contains all analyses done for the manuscript “Contrasting resistance and resilience of the coupled oxic and anoxic components of an experimental microbial community”. The first part deals with the oxygen data, the second part deals with the microbial 16S rRNA amplicons (PacBio full length sequencing).

Short summary of the project:

-8 micro-ecosystems, incubated at 24 °C for 8 days under light:dark cycles of 16:8h -afterwards, 4 micro-ecosystems (stressor treatment) were incubated in darkness for 7 days -afterwards, the stressor-treated columns were incubated again as the controls -Oxygen was measured at the top and bottom of the liquid part continuously every 5 min -samples were taken at day 8 (prior-stressor treatment), day 15 ( stressor-treatment), day 19 (short time recovery) and day 35 (long term recovery) at the height of the top and bottom oxygen sensor, respectively, for 16S rRNA full length amplicon sequencing

Aim of the study: Analysis of resistance and resilience of the oxygen concentration and the oxic/anoxic microbial community

Oxygen

data handling

oxygen_data <- read_excel("oxygen_data")

Fix the date and times

oxygen_data <- oxygen_data %>%
  mutate(Date_time = mdy_hms(paste0(Date, Time))) %>%
  select(-Date, -Time)

Calculate hourly average

hourly_value <- oxygen_data %>%
  mutate(Date = date(Date_time),
         Hour = hour(Date_time)) %>%
  group_by(Sensor_Name, Date, Hour) %>%
  dplyr::summarise(mean_Value = mean(Value)) %>%
  mutate(Date_time = ymd_h(paste0(Date, Hour)))
## `summarise()` has grouped output by 'Sensor_Name', 'Date'. You can override using the `.groups` argument.
## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.

## Warning: All formats failed to parse. No formats found.
hourly_value <- hourly_value %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")))

Calculate the DAILY mean, minimum, and maximum of the hourly means

daily_values <- hourly_value %>%
  ungroup() %>%
  group_by(Date, Sensor_Name) %>%
  dplyr::summarise(mean_O2 = mean(mean_Value),
            min_O2 = min(mean_Value),
            max_O2 = max(mean_Value),
            amplitude_O2 = max_O2 - min_O2) %>%
  mutate(Date_time_midpoint = ymd_hms(paste0(Date, "12:00:00"))) %>%
  pivot_longer(names_to = "Variable", values_to = "Oxygen", cols = 3:6)
## `summarise()` has grouped output by 'Date'. You can override using the `.groups` argument.

put in missing values for the days where there is no data

missing_data <- crossing(Sensor_Name = unique(daily_values$Sensor_Name),
                         Variable = unique(daily_values$Variable),
                         Date = dmy(paste(c(24:27), "12-2020")),
                         Date_time_midpoint = ymd_hms(paste0(Date, " 12:00:00")))%>%
  mutate(Oxygen = NA)

daily_values <- rbind(daily_values, missing_data) %>%
  mutate(Sensor_Name=tolower(Sensor_Name)) %>%
  separate(Sensor_Name, into= c("Treatment", "Sensorposition","Column"), 
           remove=FALSE)

daily_values <- daily_values %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")),
                                Variable2 = factor(Variable, levels = c("mean_O2", "max_O2", "min_O2", "amplitude_O2")))

plot the data

hourly values

plot_oxygen_a  <- hourly_value %>%
  mutate(Day2 = Day*24 + Hour) %>%
  ggplot(aes(x = Day2, y = mean_Value, col = Sensor_Name)) +
  geom_line()  +
  geom_rect(xmin = 195, xmax = 360,ymin = 0, ymax = 100, col = "lightblue", fill= "lightblue") +
  geom_line()  +
  geom_rect(xmin = 510, xmax = 640,ymin = 0, ymax = 100, col = "white", fill= "white") +
  scale_color_manual(values =c("Control_top_1" = "black", "Control_top_2" = "black", "Control_top_3" = "black", "Control_top_4" = "black","Disturbance_top_1" ="red","Disturbance_top_2" ="red","Disturbance_top_3" ="red","Disturbance_top_4" ="red","Control_bottom_1" = "grey", "Control_bottom_2" = "grey", "Control_bottom_3" = "grey", "Control_bottom_4" = "grey","Disturbance_bottom_1" ="grey","disturbance_bottom_2" ="grey","disturbance_bottom_3" ="grey","Disturbance_bottom_4" ="grey")) +
    theme_bw() + theme(panel.grid = element_blank()) + scale_x_continuous(breaks = seq(0, 817, by=48), labels=seq(0,35, by=2))+
  theme(legend.position = "none") +
  ylab("Hourly mean oxygen value [%]") +
  xlab("Day") +
  geom_point(data=data.frame(Day2=c(190, 360, 456, 840), mean_Value=rep(-1, 4)), colour="black", size=3, shape=4)


plot(plot_oxygen_a)

just top sensor with all 4 variables and daily mean

plot_oxygen_b<- daily_values %>%
  filter(Sensor_Name %in% c("control_top_1","control_top_2","control_top_3","control_top_4","disturbance_top_1","disturbance_top_2","disturbance_top_3","disturbance_top_4")) %>%
  ggplot(aes(x = Day, y = Oxygen, col = Sensor_Name)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_line()  +
  facet_wrap( ~ Variable2, ncol = 4) +
  scale_color_manual(values =c("control_top_1" = "black", "control_top_2" = "black", "control_top_3" = "black", "control_top_4" = "black","disturbance_top_1" ="red","disturbance_top_2" ="red","disturbance_top_3" ="red","disturbance_top_4" ="red","control_bottom_1" = "black", "control_bottom_2" = "black", "control_bottom_3" = "black", "control_bottom_4" = "black","disturbance_bottom_1" ="red","disturbance_bottom_2" ="red","disturbance_bottom_3" ="red","disturbance_bottom_4" ="red")) +
  theme_bw() + theme(panel.grid = element_blank())+
  theme(legend.position = "none") +
  geom_point(data=data.frame(Day=c(7.5, 15, 19, 35), Oxygen=rep(-7, 4)), shape=4 ,colour="black", size=2) +
  ylab("Daily mean oxygen value [%]") +
  xlab("Day")

plot(plot_oxygen_b)

Applying t-test for every day to get p.value for top sensors

analysis <- daily_values %>%
  filter(Sensor_Name != "control_top_1") %>%
  na.omit() %>%
  group_by(Date, Variable, Date_time_midpoint, Sensorposition) %>%
  do(m1 = tidy(lm(Oxygen ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
  unnest() %>%
  filter(term=="Treatmentdisturbance", Sensorposition == "top")
## Warning: `cols` is now required when using unnest().
## Please use `cols = c(m1)`
analysis <- analysis %>% mutate(Day = as.numeric(Date - ymd("2020-12-02")),
                                Variable2 = factor(Variable, levels = c("mean_O2", "max_O2", "min_O2", "amplitude_O2")))

Plot the estimates

#final

plot_estimate_oxygen<- analysis %>%
  ggplot(aes(x = Day, y = estimate)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point()  +
    geom_ribbon(aes(ymin=conf.low,ymax=conf.high, x=Day), alpha=.4)+
    geom_rect(xmin = 22, xmax = 25,
          ymin = -Inf, ymax = 100, col = "white", fill="white", alpha=1) +
  facet_wrap( ~ Variable2, ncol = 4) +
  scale_color_manual(values =c("control_top_1" = "black", "control_top_2" = "black", "control_top_3" = "black", "control_top_4" = "black","disturbance_top_1" ="red","disturbance_top_2" ="red","disturbance_top_3" ="red","disturbance_top_4" ="red","control_bottom_1" = "black", "control_bottom_2" = "black", "control_bottom_3" = "black", "control_bottom_4" = "black","disturbance_bottom_1" ="red","disturbance_bottom_2" ="red","disturbance_bottom_3" ="red","disturbance_bottom_4" ="red")) +
  theme_bw() + theme(panel.grid = element_blank())+ 
  geom_hline(yintercept=0, linetype="dashed") +
  geom_point(data=data.frame(Day=c(7.5, 15, 19, 35), estimate=rep(-25, 4)), colour="black", size=2, shape=4) +
  ylab("Estimate") +
  xlab("Day")

print(plot_estimate_oxygen)

Microbial community

Load phyloseq data, remove chloroplast sequences, just include samples from Experiment part III

ps <- readRDS("ps_oxic_anoxic.rds")
ps_new <-  subset_taxa(ps, Order != "Chloroplast")
taxa_names(ps_new) <- paste0("Seq", seq(ntaxa(ps)))

Calculation of alpha diversity index

diversity_test<- estimate_richness(ps_new, split=TRUE, measures=NULL)
## Warning in estimate_richness(ps_new, split = TRUE, measures = NULL): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.
diversity_test$ID <- rownames(diversity_test) 
names <- read_excel("barcodes_FKIII.xlsx")
## New names:
## * `` -> ...1
names1 <- as.data.frame(names)
names1$ID <- gsub("--", "..", names1$...1)
diversity_analysis <- merge(diversity_test, names1, by=c("ID"))

stderr <- function(x, na.rm=FALSE) {
  if (na.rm) x <- na.omit(x)
  sqrt(var(x)/length(x))
}

div_means <- diversity_analysis %>%
  dplyr::filter(Position %in% c("bottom", "top")) %>%
    dplyr::group_by(Position, Day, Treatment) %>%
    dplyr::summarize(Shannon_mean = mean(Shannon),
                     Shannon_SE = stderr(Shannon))
## `summarise()` has grouped output by 'Position', 'Day'. You can override using the `.groups` argument.
div_raw <- diversity_analysis

plot_a_diversity_top <- ggplot() +
  geom_rect(data=div_raw, aes(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = Inf), col = "lightblue", fill="lightblue", alpha=0.1) +
  geom_point(data = subset(div_raw, Position == "top"), aes(x=Day, y=Shannon, col=Treatment,shape=Position), size=3, alpha=.7)+
  scale_shape_manual(values=c(17,16))+
  geom_point(data=subset(div_means, Position == "top"), aes(x=Day, y=Shannon_mean, col=Treatment), shape="+", size = 10, position = "dodge", width = 0.25) +
  geom_line(data=subset(div_means, Position =="top"), aes(x=Day, y=Shannon_mean, col=Treatment), linetype="dashed") +
    scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
  theme_bw() + theme(panel.grid = element_blank())+ 
  theme(legend.position = "none")  +
  scale_x_continuous(limits=c(0,36),breaks = c(0,8,15,19,35)) +
  ggtitle("top layer community") +
  theme(plot.title = element_text(size = 11))
## Warning: Ignoring unknown parameters: width
print(plot_a_diversity_top)
## Warning: Width not defined. Set with `position_dodge(width = ?)`

plot_a_diversity_bottom <- ggplot() +
  geom_rect(data=div_raw, aes(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = Inf), col = "lightblue", fill="lightblue", alpha=0.1) +
  geom_point(data = subset(div_raw, Position == "bottom"), aes(x=Day, y=Shannon, col=Treatment,shape=Position), size=3, alpha=.7)+
  scale_shape_manual(values=c(17,16))+
  geom_point(data=subset(div_means, Position == "bottom"), aes(x=Day, y=Shannon_mean, col=Treatment), shape="+", size = 10, position = "dodge", width = 0.25) +
  geom_line(data=subset(div_means, Position =="bottom"), aes(x=Day, y=Shannon_mean, col=Treatment), linetype="dashed") +
    scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
  theme_bw() + theme(panel.grid = element_blank())+ 
  theme(legend.position = "none")  +
  scale_x_continuous(limits=c(0,36),breaks = c(0,8,15,19,35)) +
  ggtitle("bottom layer community")+
  theme(plot.title = element_text(size = 11))
## Warning: Ignoring unknown parameters: width
print(plot_a_diversity_bottom)
## Warning: Width not defined. Set with `position_dodge(width = ?)`

estimates of the alpha diversity analysis

 analysis_shannon <- div_raw %>%
  group_by(Day, Position) %>%
  do(m1 = tidy(lm(Shannon ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
  unnest(cols = c(m1)) %>%
  filter(term=="TreatmentDisturbance")

# plot the estimates
plot_estimates_shannon_top <- analysis_shannon%>%
  filter(Position=="top") %>%
  ggplot(aes(x = Day, y = estimate)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point()  +
    geom_errorbar(aes(ymin=conf.low,ymax=conf.high, x=Day),width=0.5, size=0.5, color="black")+
  theme_bw() + theme(panel.grid = element_blank())+ 
  geom_hline(yintercept=0, linetype="dashed") +
  ylab("Estimate") +
  xlab("Day")  +
  scale_x_continuous(limits=c(0,37),breaks = c(0, 8,15,19,35))
 

print(plot_estimates_shannon_top)

plot_estimates_shannon_bottom <- analysis_shannon%>%
  filter(Position=="bottom") %>%
  ggplot(aes(x = Day, y = estimate)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point()  +
    geom_errorbar(aes(ymin=conf.low,ymax=conf.high, x=Day),width=0.5, size=0.5, color="black")+
  theme_bw() + theme(panel.grid = element_blank())+ 
  geom_hline(yintercept=0, linetype="dashed") +
  ylab("Estimate") +
  xlab("Day")  +
  scale_x_continuous(limits=c(0,37),breaks = c(0, 8,15,19,35)) 

print(plot_estimates_shannon_bottom)

#Anova 
## top community
diversity_analysis_top<- diversity_analysis %>%
  filter(Position=="top")

lm_treatment <- lm(Shannon ~ Treatment*Day, data = diversity_analysis_top)
autoplot(lm_treatment)
## Warning: `arrange_()` was deprecated in dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help

anova(lm_treatment)
## Analysis of Variance Table
## 
## Response: Shannon
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## Treatment      1 0.0138  0.0138  0.0423 0.8385271    
## Day            1 6.5483  6.5483 20.0874 0.0001225 ***
## Treatment:Day  1 1.6860  1.6860  5.1719 0.0311192 *  
## Residuals     27 8.8018  0.3260                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm_treatment)
## 
## Call:
## lm(formula = Shannon ~ Treatment * Day, data = diversity_analysis_top)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5282 -0.3754  0.1258  0.4857  0.6453 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.51565    0.31188  17.685 2.25e-16 ***
## TreatmentDisturbance     -0.84603    0.44627  -1.896   0.0687 .  
## Day                      -0.06889    0.01441  -4.782 5.47e-05 ***
## TreatmentDisturbance:Day  0.04647    0.02044   2.274   0.0311 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.571 on 27 degrees of freedom
## Multiple R-squared:  0.4838, Adjusted R-squared:  0.4264 
## F-statistic: 8.434 on 3 and 27 DF,  p-value: 0.0004077
##top sensor shannon index is affected "postponed", and it does not show resilience
 Anova_treatment <- lm(Shannon ~ Treatment*as.factor(Day), data = diversity_analysis_top)
autoplot(Anova_treatment)

anova(Anova_treatment)
## Analysis of Variance Table
## 
## Response: Shannon
##                          Df  Sum Sq Mean Sq F value    Pr(>F)    
## Treatment                 1  0.0138  0.0138  0.0988  0.756165    
## as.factor(Day)            3 11.7919  3.9306 28.1259 7.133e-08 ***
## Treatment:as.factor(Day)  3  2.0299  0.6766  4.8418  0.009353 ** 
## Residuals                23  3.2143  0.1398                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(Anova_treatment)
## 
## Call:
## lm(formula = Shannon ~ Treatment * as.factor(Day), data = diversity_analysis_top)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68415 -0.19311  0.01744  0.30636  0.44183 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             5.3096     0.1869  28.406  < 2e-16 ***
## TreatmentDisturbance                   -0.4739     0.2643  -1.793  0.08616 .  
## as.factor(Day)15                       -0.6177     0.2643  -2.337  0.02852 *  
## as.factor(Day)19                       -1.9468     0.2643  -7.365 1.72e-07 ***
## as.factor(Day)35                       -1.9154     0.2643  -7.246 2.25e-07 ***
## TreatmentDisturbance:as.factor(Day)15   0.0423     0.3891   0.109  0.91436    
## TreatmentDisturbance:as.factor(Day)19   0.8408     0.3738   2.249  0.03438 *  
## TreatmentDisturbance:as.factor(Day)35   1.1886     0.3738   3.179  0.00418 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3738 on 23 degrees of freedom
## Multiple R-squared:  0.8115, Adjusted R-squared:  0.7541 
## F-statistic: 14.14 on 7 and 23 DF,  p-value: 5.529e-07
AIC(lm_treatment, Anova_treatment)
##                 df      AIC
## lm_treatment     5 58.94428
## Anova_treatment  9 35.71636
#bottom community
diversity_analysis_bottom<- diversity_analysis %>%
  filter(Position=="bottom")

##bottom sensor shannon index is affected directly, and it does show resilience
lm_treatment <- lm(Shannon ~ Treatment*Day, data = diversity_analysis_bottom)
autoplot(lm_treatment)

anova(lm_treatment)
## Analysis of Variance Table
## 
## Response: Shannon
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## Treatment      1  6.0279  6.0279  7.9775  0.008631 ** 
## Day            1 22.0637 22.0637 29.1999 9.209e-06 ***
## Treatment:Day  1  0.2105  0.2105  0.2786  0.601757    
## Residuals     28 21.1570  0.7556                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm_treatment)
## 
## Call:
## lm(formula = Shannon ~ Treatment * Day, data = diversity_analysis_bottom)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7771 -0.6599  0.0257  0.7631  1.4147 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               3.85624    0.47482   8.121 7.67e-09 ***
## TreatmentDisturbance      1.18319    0.67150   1.762  0.08898 .  
## Day                      -0.07561    0.02193  -3.448  0.00181 ** 
## TreatmentDisturbance:Day -0.01637    0.03102  -0.528  0.60176    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8693 on 28 degrees of freedom
## Multiple R-squared:  0.5722, Adjusted R-squared:  0.5264 
## F-statistic: 12.49 on 3 and 28 DF,  p-value: 2.307e-05
#bottom sensor shannon index is affected directly, and it does show resilience
Anova_treatment <- lm(Shannon ~ Treatment*as.factor(Day), data = diversity_analysis_bottom)
autoplot(Anova_treatment)

anova(Anova_treatment)
## Analysis of Variance Table
## 
## Response: Shannon
##                          Df  Sum Sq Mean Sq F value    Pr(>F)    
## Treatment                 1  6.0279  6.0279  12.050  0.001978 ** 
## as.factor(Day)            3 25.2738  8.4246  16.841 4.204e-06 ***
## Treatment:as.factor(Day)  3  6.1515  2.0505   4.099  0.017537 *  
## Residuals                24 12.0059  0.5002                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(Anova_treatment)
## 
## Call:
## lm(formula = Shannon ~ Treatment * as.factor(Day), data = diversity_analysis_bottom)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.22400 -0.31305  0.00521  0.39579  1.33089 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            3.39590    0.35364   9.603 1.08e-09 ***
## TreatmentDisturbance                  -0.06541    0.50012  -0.131 0.897026    
## as.factor(Day)15                      -0.73628    0.50012  -1.472 0.153962    
## as.factor(Day)19                      -1.14224    0.50012  -2.284 0.031520 *  
## as.factor(Day)35                      -2.10233    0.50012  -4.204 0.000315 ***
## TreatmentDisturbance:as.factor(Day)15  1.83580    0.70728   2.596 0.015862 *  
## TreatmentDisturbance:as.factor(Day)19  1.78258    0.70728   2.520 0.018783 *  
## TreatmentDisturbance:as.factor(Day)35  0.11542    0.70728   0.163 0.871739    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7073 on 24 degrees of freedom
## Multiple R-squared:  0.7573, Adjusted R-squared:  0.6865 
## F-statistic:  10.7 on 7 and 24 DF,  p-value: 4.674e-06
AIC(lm_treatment, Anova_treatment)
##                 df      AIC
## lm_treatment     5 87.57164
## Anova_treatment  9 77.44136

Calculation rel abundance, just sequences that appear > 0.1 %, making dataframes for each column position

ps_rel <- transform_sample_counts(ps_new, function(x) x / sum(x) )

relab_threshold <- 0.001

ps_relab <- filter_taxa(ps_rel, function(x) !(sum(x < relab_threshold) == length(x)), TRUE)
ntaxa(ps_new)
## [1] 9028
ntaxa(ps_relab)
## [1] 1856
ps_relative <- transform_sample_counts(ps_relab, function(x)  x / sum(x))
df_rel <- psmelt(ps_relative)
df_rel_top <- df_rel %>%
  filter(Position == "top")
df_rel_bottom <- df_rel %>%
  filter(Position == "bottom")

NMDS

NMDS Plot

mds_whole <- ps_relative@otu_table %>%
  as.data.frame() %>%
  metaMDS(., 
        distance = "bray",
        k = 2, ## number of dimensions to reduce to
        try = 200, ## number of random starts to try
        autotransform = FALSE ## best not to use
)
## Run 0 stress 0.1530756 
## Run 1 stress 0.1530827 
## ... Procrustes: rmse 0.00148772  max resid 0.0103001 
## Run 2 stress 0.1530831 
## ... Procrustes: rmse 0.001531195  max resid 0.01032903 
## Run 3 stress 0.1571977 
## Run 4 stress 0.236323 
## Run 5 stress 0.1571977 
## Run 6 stress 0.1571977 
## Run 7 stress 0.1530827 
## ... Procrustes: rmse 0.001485717  max resid 0.01028305 
## Run 8 stress 0.1530829 
## ... Procrustes: rmse 0.001513705  max resid 0.01031943 
## Run 9 stress 0.1530827 
## ... Procrustes: rmse 0.001483084  max resid 0.01026715 
## Run 10 stress 0.1530827 
## ... Procrustes: rmse 0.00148431  max resid 0.01027275 
## Run 11 stress 0.2078723 
## Run 12 stress 0.1530827 
## ... Procrustes: rmse 0.001483782  max resid 0.01027015 
## Run 13 stress 0.153083 
## ... Procrustes: rmse 0.001527468  max resid 0.01032705 
## Run 14 stress 0.153083 
## ... Procrustes: rmse 0.001522597  max resid 0.0103238 
## Run 15 stress 0.1530827 
## ... Procrustes: rmse 0.001484558  max resid 0.01027117 
## Run 16 stress 0.1530827 
## ... Procrustes: rmse 0.001486136  max resid 0.01027742 
## Run 17 stress 0.1571977 
## Run 18 stress 0.1530827 
## ... Procrustes: rmse 0.001486199  max resid 0.01028388 
## Run 19 stress 0.1571977 
## Run 20 stress 0.1571814 
## Run 21 stress 0.1571977 
## Run 22 stress 0.2275183 
## Run 23 stress 0.1530827 
## ... Procrustes: rmse 0.00148717  max resid 0.01029209 
## Run 24 stress 0.1530827 
## ... Procrustes: rmse 0.00148708  max resid 0.0102902 
## Run 25 stress 0.1835705 
## Run 26 stress 0.1530829 
## ... Procrustes: rmse 0.001522374  max resid 0.01032379 
## Run 27 stress 0.1530756 
## ... Procrustes: rmse 2.585429e-05  max resid 0.0001694419 
## ... Similar to previous best
## Run 28 stress 0.1530827 
## ... Procrustes: rmse 0.001485381  max resid 0.01027804 
## Run 29 stress 0.1863566 
## Run 30 stress 0.1571814 
## Run 31 stress 0.1530756 
## ... New best solution
## ... Procrustes: rmse 1.911943e-05  max resid 0.0001243164 
## ... Similar to previous best
## Run 32 stress 0.1530756 
## ... Procrustes: rmse 4.391557e-05  max resid 0.0002895892 
## ... Similar to previous best
## Run 33 stress 0.2462748 
## Run 34 stress 0.2452603 
## Run 35 stress 0.1530756 
## ... Procrustes: rmse 2.010439e-05  max resid 0.0001319844 
## ... Similar to previous best
## Run 36 stress 0.1530827 
## ... Procrustes: rmse 0.00148587  max resid 0.01028955 
## Run 37 stress 0.1530756 
## ... Procrustes: rmse 1.716225e-05  max resid 0.0001121669 
## ... Similar to previous best
## Run 38 stress 0.1571977 
## Run 39 stress 0.227854 
## Run 40 stress 0.1571977 
## Run 41 stress 0.1530827 
## ... Procrustes: rmse 0.001486417  max resid 0.01029439 
## Run 42 stress 0.1530756 
## ... Procrustes: rmse 1.800349e-05  max resid 0.0001184586 
## ... Similar to previous best
## Run 43 stress 0.2437621 
## Run 44 stress 0.153083 
## ... Procrustes: rmse 0.001520742  max resid 0.01031471 
## Run 45 stress 0.1530756 
## ... Procrustes: rmse 2.018693e-05  max resid 0.0001329001 
## ... Similar to previous best
## Run 46 stress 0.1571977 
## Run 47 stress 0.1530756 
## ... Procrustes: rmse 2.90842e-06  max resid 1.408669e-05 
## ... Similar to previous best
## Run 48 stress 0.1530827 
## ... Procrustes: rmse 0.001482603  max resid 0.01026852 
## Run 49 stress 0.1530827 
## ... Procrustes: rmse 0.001485844  max resid 0.01028957 
## Run 50 stress 0.1530827 
## ... Procrustes: rmse 0.00148322  max resid 0.01027252 
## Run 51 stress 0.1530827 
## ... Procrustes: rmse 0.001485107  max resid 0.01028072 
## Run 52 stress 0.1530829 
## ... Procrustes: rmse 0.001515176  max resid 0.01032225 
## Run 53 stress 0.1530831 
## ... Procrustes: rmse 0.001532307  max resid 0.0103358 
## Run 54 stress 0.153083 
## ... Procrustes: rmse 0.001529385  max resid 0.01033284 
## Run 55 stress 0.1530827 
## ... Procrustes: rmse 0.001482715  max resid 0.01027004 
## Run 56 stress 0.183571 
## Run 57 stress 0.1571977 
## Run 58 stress 0.153083 
## ... Procrustes: rmse 0.001501185  max resid 0.01030609 
## Run 59 stress 0.1530756 
## ... Procrustes: rmse 2.806648e-06  max resid 1.687471e-05 
## ... Similar to previous best
## Run 60 stress 0.1530827 
## ... Procrustes: rmse 0.001486085  max resid 0.01028057 
## Run 61 stress 0.1530831 
## ... Procrustes: rmse 0.00154512  max resid 0.01034472 
## Run 62 stress 0.1571977 
## Run 63 stress 0.1530756 
## ... New best solution
## ... Procrustes: rmse 9.745351e-06  max resid 6.388499e-05 
## ... Similar to previous best
## Run 64 stress 0.1530827 
## ... Procrustes: rmse 0.001483923  max resid 0.01027472 
## Run 65 stress 0.2383096 
## Run 66 stress 0.1530757 
## ... Procrustes: rmse 0.0001074537  max resid 0.0007097788 
## ... Similar to previous best
## Run 67 stress 0.2078723 
## Run 68 stress 0.153083 
## ... Procrustes: rmse 0.001533951  max resid 0.01033842 
## Run 69 stress 0.1530828 
## ... Procrustes: rmse 0.001502232  max resid 0.01030918 
## Run 70 stress 0.1530827 
## ... Procrustes: rmse 0.00148633  max resid 0.01028869 
## Run 71 stress 0.1530827 
## ... Procrustes: rmse 0.001487739  max resid 0.01028381 
## Run 72 stress 0.1571977 
## Run 73 stress 0.1530756 
## ... Procrustes: rmse 1.115265e-05  max resid 7.306532e-05 
## ... Similar to previous best
## Run 74 stress 0.1530832 
## ... Procrustes: rmse 0.001553023  max resid 0.01035205 
## Run 75 stress 0.1530827 
## ... Procrustes: rmse 0.001481748  max resid 0.0102606 
## Run 76 stress 0.1530828 
## ... Procrustes: rmse 0.001498745  max resid 0.01030342 
## Run 77 stress 0.1571977 
## Run 78 stress 0.1530827 
## ... Procrustes: rmse 0.001488154  max resid 0.01029276 
## Run 79 stress 0.1530827 
## ... Procrustes: rmse 0.001492443  max resid 0.01029632 
## Run 80 stress 0.1530827 
## ... Procrustes: rmse 0.001488339  max resid 0.01030626 
## Run 81 stress 0.1530756 
## ... Procrustes: rmse 1.113055e-05  max resid 7.311672e-05 
## ... Similar to previous best
## Run 82 stress 0.1530829 
## ... Procrustes: rmse 0.001509633  max resid 0.0103157 
## Run 83 stress 0.1530756 
## ... Procrustes: rmse 7.169728e-06  max resid 4.657897e-05 
## ... Similar to previous best
## Run 84 stress 0.1530827 
## ... Procrustes: rmse 0.00148912  max resid 0.01028799 
## Run 85 stress 0.226658 
## Run 86 stress 0.1530828 
## ... Procrustes: rmse 0.001508742  max resid 0.01031494 
## Run 87 stress 0.1530827 
## ... Procrustes: rmse 0.00148184  max resid 0.01026222 
## Run 88 stress 0.2288104 
## Run 89 stress 0.1530831 
## ... Procrustes: rmse 0.001544824  max resid 0.01035264 
## Run 90 stress 0.1530831 
## ... Procrustes: rmse 0.001545692  max resid 0.01035093 
## Run 91 stress 0.1530831 
## ... Procrustes: rmse 0.001537873  max resid 0.0103442 
## Run 92 stress 0.1530756 
## ... Procrustes: rmse 2.447473e-06  max resid 1.156627e-05 
## ... Similar to previous best
## Run 93 stress 0.1530827 
## ... Procrustes: rmse 0.001483834  max resid 0.01027246 
## Run 94 stress 0.1530828 
## ... Procrustes: rmse 0.00150334  max resid 0.01031119 
## Run 95 stress 0.2383482 
## Run 96 stress 0.1530829 
## ... Procrustes: rmse 0.00151533  max resid 0.01032149 
## Run 97 stress 0.153083 
## ... Procrustes: rmse 0.00153258  max resid 0.01033407 
## Run 98 stress 0.1530827 
## ... Procrustes: rmse 0.001490857  max resid 0.01029209 
## Run 99 stress 0.1530829 
## ... Procrustes: rmse 0.001515011  max resid 0.01032188 
## Run 100 stress 0.2044065 
## Run 101 stress 0.2284918 
## Run 102 stress 0.1530827 
## ... Procrustes: rmse 0.00148422  max resid 0.01027434 
## Run 103 stress 0.1835701 
## Run 104 stress 0.1530827 
## ... Procrustes: rmse 0.001487628  max resid 0.01029498 
## Run 105 stress 0.1530829 
## ... Procrustes: rmse 0.001519059  max resid 0.01032514 
## Run 106 stress 0.1530827 
## ... Procrustes: rmse 0.001482414  max resid 0.010266 
## Run 107 stress 0.1571977 
## Run 108 stress 0.2233197 
## Run 109 stress 0.153083 
## ... Procrustes: rmse 0.00153533  max resid 0.01034335 
## Run 110 stress 0.1530828 
## ... Procrustes: rmse 0.001504003  max resid 0.01031102 
## Run 111 stress 0.1571977 
## Run 112 stress 0.2198539 
## Run 113 stress 0.1530827 
## ... Procrustes: rmse 0.001487103  max resid 0.01029527 
## Run 114 stress 0.1530829 
## ... Procrustes: rmse 0.001510092  max resid 0.01031824 
## Run 115 stress 0.1530756 
## ... Procrustes: rmse 5.588646e-06  max resid 3.533489e-05 
## ... Similar to previous best
## Run 116 stress 0.1530757 
## ... Procrustes: rmse 8.080733e-05  max resid 0.0005236894 
## ... Similar to previous best
## Run 117 stress 0.1530756 
## ... Procrustes: rmse 1.849313e-06  max resid 1.105e-05 
## ... Similar to previous best
## Run 118 stress 0.1530756 
## ... Procrustes: rmse 2.52858e-06  max resid 1.035051e-05 
## ... Similar to previous best
## Run 119 stress 0.1571977 
## Run 120 stress 0.1530828 
## ... Procrustes: rmse 0.001505058  max resid 0.01031287 
## Run 121 stress 0.1571977 
## Run 122 stress 0.236946 
## Run 123 stress 0.2201603 
## Run 124 stress 0.1571814 
## Run 125 stress 0.1530831 
## ... Procrustes: rmse 0.001543113  max resid 0.0103587 
## Run 126 stress 0.1530756 
## ... Procrustes: rmse 8.588398e-06  max resid 5.180093e-05 
## ... Similar to previous best
## Run 127 stress 0.1835683 
## Run 128 stress 0.1571977 
## Run 129 stress 0.1530756 
## ... Procrustes: rmse 1.594459e-06  max resid 4.629072e-06 
## ... Similar to previous best
## Run 130 stress 0.1530831 
## ... Procrustes: rmse 0.001539571  max resid 0.01036053 
## Run 131 stress 0.1530831 
## ... Procrustes: rmse 0.001537212  max resid 0.01033847 
## Run 132 stress 0.1530827 
## ... Procrustes: rmse 0.001483814  max resid 0.01027562 
## Run 133 stress 0.2044065 
## Run 134 stress 0.1530756 
## ... Procrustes: rmse 1.198952e-05  max resid 7.876109e-05 
## ... Similar to previous best
## Run 135 stress 0.2078903 
## Run 136 stress 0.1864996 
## Run 137 stress 0.1571977 
## Run 138 stress 0.1530827 
## ... Procrustes: rmse 0.001483878  max resid 0.01027624 
## Run 139 stress 0.1530756 
## ... Procrustes: rmse 1.102587e-05  max resid 7.150153e-05 
## ... Similar to previous best
## Run 140 stress 0.2418881 
## Run 141 stress 0.2149961 
## Run 142 stress 0.2247954 
## Run 143 stress 0.1571814 
## Run 144 stress 0.1571814 
## Run 145 stress 0.1530827 
## ... Procrustes: rmse 0.00148572  max resid 0.01028685 
## Run 146 stress 0.1571977 
## Run 147 stress 0.1530827 
## ... Procrustes: rmse 0.00148448  max resid 0.01027579 
## Run 148 stress 0.1530829 
## ... Procrustes: rmse 0.001512976  max resid 0.01031603 
## Run 149 stress 0.1835696 
## Run 150 stress 0.1571977 
## Run 151 stress 0.2202608 
## Run 152 stress 0.1530756 
## ... Procrustes: rmse 8.212996e-06  max resid 5.359464e-05 
## ... Similar to previous best
## Run 153 stress 0.1530756 
## ... Procrustes: rmse 7.548606e-06  max resid 4.950594e-05 
## ... Similar to previous best
## Run 154 stress 0.1530827 
## ... Procrustes: rmse 0.001486845  max resid 0.01029449 
## Run 155 stress 0.1530827 
## ... Procrustes: rmse 0.001480886  max resid 0.01024511 
## Run 156 stress 0.1530829 
## ... Procrustes: rmse 0.001516893  max resid 0.010325 
## Run 157 stress 0.1530827 
## ... Procrustes: rmse 0.001479907  max resid 0.01024182 
## Run 158 stress 0.2203492 
## Run 159 stress 0.1530827 
## ... Procrustes: rmse 0.001487245  max resid 0.01029409 
## Run 160 stress 0.1530827 
## ... Procrustes: rmse 0.001485919  max resid 0.01028685 
## Run 161 stress 0.1530757 
## ... Procrustes: rmse 0.0001006575  max resid 0.0006602985 
## ... Similar to previous best
## Run 162 stress 0.1530827 
## ... Procrustes: rmse 0.001487802  max resid 0.01029371 
## Run 163 stress 0.1530829 
## ... Procrustes: rmse 0.00150922  max resid 0.01030585 
## Run 164 stress 0.1530828 
## ... Procrustes: rmse 0.001495956  max resid 0.01030119 
## Run 165 stress 0.2348738 
## Run 166 stress 0.153083 
## ... Procrustes: rmse 0.001530306  max resid 0.01033334 
## Run 167 stress 0.1571977 
## Run 168 stress 0.1571977 
## Run 169 stress 0.1530831 
## ... Procrustes: rmse 0.001540548  max resid 0.01034203 
## Run 170 stress 0.1530827 
## ... Procrustes: rmse 0.001483741  max resid 0.01027234 
## Run 171 stress 0.1530756 
## ... Procrustes: rmse 4.863105e-06  max resid 3.150449e-05 
## ... Similar to previous best
## Run 172 stress 0.1571977 
## Run 173 stress 0.1530756 
## ... Procrustes: rmse 1.542525e-06  max resid 7.538724e-06 
## ... Similar to previous best
## Run 174 stress 0.2044065 
## Run 175 stress 0.153083 
## ... Procrustes: rmse 0.001528137  max resid 0.01033225 
## Run 176 stress 0.1530829 
## ... Procrustes: rmse 0.001509756  max resid 0.01031674 
## Run 177 stress 0.1530831 
## ... Procrustes: rmse 0.001545012  max resid 0.01034572 
## Run 178 stress 0.2427732 
## Run 179 stress 0.1530827 
## ... Procrustes: rmse 0.001492616  max resid 0.01029683 
## Run 180 stress 0.1530827 
## ... Procrustes: rmse 0.001483253  max resid 0.01027038 
## Run 181 stress 0.2427561 
## Run 182 stress 0.1530831 
## ... Procrustes: rmse 0.001547099  max resid 0.01034647 
## Run 183 stress 0.2198537 
## Run 184 stress 0.1571814 
## Run 185 stress 0.1530756 
## ... Procrustes: rmse 9.468755e-06  max resid 6.206903e-05 
## ... Similar to previous best
## Run 186 stress 0.1571977 
## Run 187 stress 0.1530756 
## ... New best solution
## ... Procrustes: rmse 1.888246e-06  max resid 1.151079e-05 
## ... Similar to previous best
## Run 188 stress 0.1530756 
## ... Procrustes: rmse 1.28233e-05  max resid 8.444567e-05 
## ... Similar to previous best
## Run 189 stress 0.1530757 
## ... Procrustes: rmse 5.772538e-05  max resid 0.0003430395 
## ... Similar to previous best
## Run 190 stress 0.2277348 
## Run 191 stress 0.2450719 
## Run 192 stress 0.2220019 
## Run 193 stress 0.1530827 
## ... Procrustes: rmse 0.001486224  max resid 0.01028027 
## Run 194 stress 0.2430565 
## Run 195 stress 0.1530828 
## ... Procrustes: rmse 0.001507977  max resid 0.0103144 
## Run 196 stress 0.153083 
## ... Procrustes: rmse 0.001525161  max resid 0.01033065 
## Run 197 stress 0.2044065 
## Run 198 stress 0.1530827 
## ... Procrustes: rmse 0.001486139  max resid 0.01028887 
## Run 199 stress 0.2373475 
## Run 200 stress 0.1571977 
## *** Solution reached
## 0.16


mds_whole_res <- ps_relative @sam_data %>%
  as.tibble() %>%
  select(Treatment, Column, Day, Position,Day_name) %>%
  bind_cols(as.tibble(scores(mds_whole, display = "sites")))
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 object
nmds_plot_supplements <- ggplot(mds_whole_res, aes(x = NMDS1, y = NMDS2)) +
  geom_point(aes(shape = Position, color=Treatment),
             size = 2) +
  scale_shape_manual(values =c(17,16))+
  facet_wrap("Day", ncol=2)+
   theme_bw() + theme(panel.grid = element_blank())+ 
  scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))

print(nmds_plot_supplements)

NMDS separated for NMDS1 and NMDS2 axis

NMDS_top <- mds_whole_res %>% 
 filter(Position=="top") 

NMDS_bottom <- mds_whole_res %>% 
 filter(Position=="bottom") 


# Dynamics
plot_NMDS1_top <- NMDS_top %>%
  ggplot(aes(x = Day, y = NMDS1, col = Treatment)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point(shape=16, size = 2)+
  geom_smooth() +
  theme_bw() + theme(panel.grid = element_blank())+ 
    scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
  ggtitle("top layer community") +
  scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
  theme(legend.position = "none")+
  theme(plot.title = element_text(size = 11))

print(plot_NMDS1_top)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42

plot_NMDS2_top <- NMDS_top %>%
  ggplot(aes(x = Day, y = NMDS2, col = Treatment)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point(shape=16, size = 2) +
  geom_smooth() +
  theme_bw() + theme(panel.grid = element_blank())+ 
  ggtitle("top layer community")+ 
  scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
    scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
  theme(legend.position = "none")+
  theme(plot.title = element_text(size = 11))

print(plot_NMDS2_top)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42

plot_NMDS1_bottom <- NMDS_bottom %>%
  ggplot(aes(x = Day, y = NMDS1, col = Treatment)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point(shape=17, size = 2) +
  geom_smooth() +
  theme_bw() + theme(panel.grid = element_blank())+ 
    scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
  ggtitle("bottom layer community")+ 
  scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
  theme(legend.position = "none")+
  theme(plot.title = element_text(size = 11))

print(plot_NMDS1_bottom)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42

plot_NMDS2_bottom<- NMDS_bottom %>%
  ggplot(aes(x = Day, y = NMDS2, col = Treatment)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point(shape=17, size = 2) +
  geom_smooth() +
  theme_bw() + theme(panel.grid = element_blank())+ 
  scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
    scale_colour_manual(values = c("Control" = "black", "Disturbance" = "red"))+
  ggtitle("bottom layer community")+ 
  theme(legend.position = "none")+
  theme(plot.title = element_text(size = 11))

print(plot_NMDS2_bottom)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42

estimate analysis for NMDS data

analysis_NMDS1_top <- NMDS_top %>%
  group_by(Day, Day_name) %>%
  do(m1 = tidy(lm(NMDS1 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
  unnest(cols = c(m1)) %>%
  filter(term=="TreatmentDisturbance")

plot_estimate_NMDS1_top <- analysis_NMDS1_top %>%
  ggplot(aes(x = Day, y = estimate)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point()  +
  geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+ 
  theme_bw() + theme(panel.grid = element_blank())+ 
  geom_hline(yintercept=0, linetype="dashed") +
  scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
  ylab("Estimate") +
  xlab("Day")

print(plot_estimate_NMDS1_top)

analysis_NMDS2_top <- NMDS_top %>%
  group_by(Day, Day_name) %>%
  do(m1 = tidy(lm(NMDS2 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
  unnest(cols = c(m1)) %>%
  filter(term=="TreatmentDisturbance")

plot_estimate_NMDS2_top <- analysis_NMDS2_top %>%
  ggplot(aes(x = Day, y = estimate)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point()  +
  geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
  theme_bw() + theme(panel.grid = element_blank())+ 
  scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
  geom_hline(yintercept=0, linetype="dashed") +
  ylab("Estimate") +
  xlab("Day")

print(plot_estimate_NMDS2_top)

analysis_NMDS1_bottom <- NMDS_bottom %>%
  group_by(Day, Day_name) %>%
  do(m1 = tidy(lm(NMDS1 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
  unnest(cols = c(m1)) %>%
  filter(term=="TreatmentDisturbance")

plot_estimate_NMDS1_bottom <- analysis_NMDS1_bottom%>%
  ggplot(aes(x = Day, y = estimate)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point()  +
  geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
  theme_bw() + theme(panel.grid = element_blank())+ 
  geom_hline(yintercept=0, linetype="dashed")+
  scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
  ylab("Estimate") +
  xlab("Day")

print(plot_estimate_NMDS1_bottom)

analysis_NMDS2_bottom <- NMDS_bottom %>%
  group_by(Day, Day_name) %>%
  do(m1 = tidy(lm(NMDS2 ~ Treatment, data=.), conf.int=T, conf.level = 0.95)) %>%
  unnest(cols = c(m1)) %>%
  filter(term=="TreatmentDisturbance")

plot_estimate_NMDS2_bottom <- analysis_NMDS2_bottom%>%
  ggplot(aes(x = Day, y = estimate)) +
  geom_rect(xmin = 8, xmax = 15,
          ymin = -Inf, ymax = 100, col = "lightblue", fill="lightblue", alpha=.1) +
  geom_point()  +
  geom_errorbar(aes(x=Day, ymin=conf.low, ymax=conf.high), width=0.5, size=0.5, color="black")+
  theme_bw() + theme(panel.grid = element_blank())+ 
  geom_hline(yintercept=0, linetype="dashed") +
  scale_x_continuous(limits=c(0,36),breaks = c(0, 8,15,19,35))+
  ylab("Estimate") +
  xlab("Day")

print(plot_estimate_NMDS2_bottom)

final plot of figures combined

p1 <- plot_oxygen_a
p2 <- plot_oxygen_b
p3 <- plot_a_diversity_top
p4<-  plot_a_diversity_bottom
p5 <- plot_estimates_shannon_top
p6 <- plot_estimates_shannon_bottom
p7 <- plot_NMDS1_bottom
p8 <- plot_NMDS2_bottom
p9 <- plot_NMDS1_top
p10 <- plot_NMDS2_top

p11 <- plot_estimate_oxygen
p12<- plot_estimates_shannon_top
p13<- plot_estimates_shannon_bottom
p14 <- plot_estimate_NMDS1_bottom
p15 <- plot_estimate_NMDS2_bottom
p16 <- plot_estimate_NMDS1_top
p17 <- plot_estimate_NMDS2_top


#Figure 1
plot_1 <- p1/
  p2/
  p11 

plot_1 +
  plot_layout(heights=c(2,1,1)) + plot_annotation(tag_levels = "a")

#Figure 2
plot <- (p4|p13)/
  (p3|p12)/
  (p7|p14)/
  (p8|p15)/
  (p9 |p16)/
  (p10 |p17)
  

plot + plot_annotation(tag_levels = "a")
## Warning: Width not defined. Set with `position_dodge(width = ?)`

## Warning: Width not defined. Set with `position_dodge(width = ?)`
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 5.4432e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 5.4432e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 7.865
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 11.135
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.8973e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 405.42
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 7.865
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 11.135
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.8973e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 405.42

#Supplementary Information

Visualization of rel.abundance for supplementary information

df_class<- df_rel %>%
  group_by(Class, Treatment,Column,Day, Position) %>%
  dplyr::summarize(Abundance=sum(Abundance))
## `summarise()` has grouped output by 'Class', 'Treatment', 'Column', 'Day'. You can override using the `.groups` argument.
index <- which(df_class$Abundance>=0.05)
class_to_keep <- unique(df_class[index,"Class"])
class_to_keep <- unname(unlist(class_to_keep))

 

df_class$Class_filter <- ifelse(df_class$Class %in% class_to_keep, df_class$Class,"Zother")

   cols <-c("Actinobacteria"="bisque1", "Alphaproteobacteria" = "yellow", "Babeliae"="lightcyan1","Bacilli" = "green4",  "Bacteroidia" ="maroon1", "Campylobacteria" ="royalblue1","Zother"="black","Chlorobia" ="thistle1", "Clostridia"= "firebrick1", "Cyanobacteriia" = "magenta4", "Fimbriimonadia" = "grey","Gammaproteobacteria"="red", "Microgenomatia"="grey38","Oligoflexia"="snow1", "Phycisphaerae"="cornflowerblue", "Planctomycetes"="lightskyblue1","Spirochaetia"="lavender", "Vampirivibronia" = "lightsalmon", "Verrucomicrobiae" = "slateblue3")
   
  supplement_plot1 <-df_class%>%
  ggplot(aes_string(x = "Column", y = "Abundance", fill="Class_filter" )) +
     geom_bar(stat = "identity", position = "stack", col="black") +
       facet_wrap(Position~Day, ncol = 4)  +
       theme(axis.text.x = element_text(angle = -45, hjust = 0, vjust=0.5, size=9))+
    ylab("Relative abundance") +
  scale_fill_manual("legend", values = cols)
  
  print(supplement_plot1)

Correlation last day, shannon vs.oxygen amplitude

new <- div_raw %>%
filter(Day=="35", Position=="top")

new2 <- daily_values %>%
  filter(Day=="35", Variable=="amplitude_O2", Sensorposition=="top") %>% mutate(Column2 = paste0(str_to_title(Treatment),"_", Column))
  

new3 <- merge(new, new2, by.x=c("Column"), by.y=c("Column2")) 



ggplot(new3) + geom_point(aes(x=log(Shannon),y=Oxygen, colour=Treatment.x),size=3, cex.lab=10) +
    theme_bw() + theme(panel.grid = element_blank()) +
  scale_color_manual(values =c("Control" = "black", "Disturbance" = "red")) +
  stat_smooth(aes(x=log(Shannon),y=Oxygen), method="lm") +
  ylab("Oxygen (Amplitude)") +
  xlab("log(Shannon-index)") +
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(size=16),axis.text.y =element_text(size=20)) +
  theme(axis.title = element_text(size=16))
## Warning: Ignoring unknown parameters: cex.lab
## `geom_smooth()` using formula 'y ~ x'

cor.test(log(new3$Shannon),new3$Oxygen)
## 
##  Pearson's product-moment correlation
## 
## data:  log(new3$Shannon) and new3$Oxygen
## t = -5.1134, df = 6, p-value = 0.002192
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9822793 -0.5410140
## sample estimates:
##        cor 
## -0.9018627
#without column 3

new3_no_control_3 <- new3 %>%
  filter(Column != "Control_3")

ggplot(new3_no_control_3) + geom_point(aes(x=log(Shannon),y=Oxygen, colour=Treatment.x),size=4, cex.lab=10) +
    theme_bw() + theme(panel.grid = element_blank()) +
  scale_color_manual(values =c("Control" = "black", "Disturbance" = "red")) +
  stat_smooth(aes(x=log(Shannon),y=Oxygen), method="lm") +
  ylab("Oxygen (Amplitude)") +
  xlab("log(Shannon-index)") +
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(size=16),axis.text.y =element_text(size=20)) +
  theme(axis.title = element_text(size=16))
## Warning: Ignoring unknown parameters: cex.lab
## `geom_smooth()` using formula 'y ~ x'

cor.test(log(new3_no_control_3$Shannon),new3_no_control_3$Oxygen)
## 
##  Pearson's product-moment correlation
## 
## data:  log(new3_no_control_3$Shannon) and new3_no_control_3$Oxygen
## t = -2.6205, df = 5, p-value = 0.04707
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.96242875 -0.01788358
## sample estimates:
##        cor 
## -0.7606971
nmds <- mds_whole_res %>%
  filter(Day=="35", Position=="top")

nmds2 <- merge(nmds, new2, by.x=c("Column"), by.y=c("Column2"))


ggplot(nmds2, aes(x=NMDS1,y=Oxygen, colour=Treatment.x)) + geom_point()

cor.test(nmds2$NMDS1,nmds2$Oxygen)
## 
##  Pearson's product-moment correlation
## 
## data:  nmds2$NMDS1 and nmds2$Oxygen
## t = 0.45373, df = 6, p-value = 0.666
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5994779  0.7859367
## sample estimates:
##       cor 
## 0.1821357
ggplot(nmds2, aes(x=NMDS2,y=Oxygen, colour=Treatment.x)) + geom_point()

cor.test(nmds2$NMDS2,nmds2$Oxygen)
## 
##  Pearson's product-moment correlation
## 
## data:  nmds2$NMDS2 and nmds2$Oxygen
## t = 1.1638, df = 6, p-value = 0.2887
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3949661  0.8705568
## sample estimates:
##       cor 
## 0.4291481